aa[rep] = test_roc$auc
roc[[rep]] = test_roc
}
mean(aa)
sd(aa)
plot(roc[[1]], col = 1, main = "ROC curves on 50 replications", lwd=0.5)
for(rep in 2:50){
plot(roc[[rep]], col=1, add=T, lwd=0.5)
}
roc[[1]]
roc[[1]]$sensitivities
roc[[1]]$specificities
res = as.data.frame(0, ncol=2)
res
res = c(roc[[1]]$specificities, roc[[1]]$sensitivities)
res
for(i in 1:50){
print(roc[[i]]$specificities)
}
aa =c()
roc = list()
for(rep in 1:50){
set.seed(rep+1234)
inTraining <- createDataPartition(dat$shift_or_not, p = .80, list = FALSE)
training <- dat[inTraining,]
testing  <- dat[-inTraining,]
training = na.omit(training)
testing = na.omit(testing)
final.mod = glm(shift_or_not ~ fh_polity2 + atop_number + bicc_milper,
data=training,
family = binomial(link="logit"))
test_prob = predict(final.mod , newdata = testing, type = "response")
test_roc = roc(testing$shift_or_not ~ test_prob, plot = FALSE, print.auc = FALSE)
# summary(final.mod)
aa[rep] = test_roc$auc
roc[[rep]] = test_roc
}
mean(aa)
sd(aa)
which.max(aa)
which.max(aa)
which.min(aa)
which(aa= median(aa))
which(aa== median(aa))
median(aa)
which(aa) ==median(aa)
aaa
aa
which.max(aa)
which.min(aa)
median(aa) #9
# median(aa) #9
plot(roc[[9]], col="blue")
which.max(aa) # 45
which.min(aa) # 30
# median(aa) #9
plot(roc[[9]], col="blue")
plot(roc[[45]], col="green", add=TRUE)
which.max(aa) # 45
which.min(aa) # 30
# median(aa) #9
plot(roc[[9]], col="blue")
plot(roc[[45]], col="green", add=TRUE)
plot(roc[[30]], col="red", add=TRUE)
which.max(aa) # 45
which.min(aa) # 30
# median(aa) #9
plot(roc[[9]], col="blue")
plot(roc[[45]], col="darkgreen", add=TRUE)
plot(roc[[30]], col="red", add=TRUE)
quantile(aa)
# median(aa) #9
q = quantile(aa)
a[2]
q[2]
which(aa == q[2])
aa
sort(aa)
order(aa)
order(aa)
order(aa) == 1
which(order(aa) == 1)
50/4
aa
order(aa)
which(order(aa) == 13)
aa[which(order(aa) == 13)]
aa[which(order(aa) == 50)]
aa
order(aa)
which(order(aa) == 50)
aa[which(order(aa) == 50)]
# median(aa) #9
q = quantile(aa)
q
summary(aa)
which(aa == 0.8137)
which(aa < 0.8137)
which.min(aa[hich(aa > 0.8137)])
which.min(aa[which(aa > 0.8137)])
aa[22]
aa[which(aa > 0.8137)]
aa
which.max(aa) # 45
which.min(aa) # 30
# median(aa) #9
q = quantile(aa)
which(aa == q[2])
plot(roc[[9]], col="red")
plot(roc[[45]], col="darkgreen", add=TRUE)
plot(roc[[30]], col="red", add=TRUE)
q = quantile(aa)
q
which.max(aa) # 45
which.min(aa) # 30
# median(aa) #9
q = quantile(aa)
which(aa == q[2])
plot(roc[[9]], col="blue")
plot(roc[[28]], col="darkgreen", add=TRUE)
plot(roc[[43]], col="red", add=TRUE)
which.max(aa) # 45
which.min(aa) # 30
# median(aa) #9
q = quantile(aa)
which(aa == q[2])
plot(roc[[9]], col="blue")
plot(roc[[28]], col="black", add=TRUE)
plot(roc[[43]], col="red", add=TRUE)
which.max(aa) # 45
which.min(aa) # 30
# median(aa) #9
q = quantile(aa)
which(aa == q[2])
plot(roc[[9]], col="blue")
plot(roc[[28]], col="black", add=TRUE, lty =2)
plot(roc[[43]], col="red", add=TRUE, lty = 2)
plot(roc[[9]], col="blue")
plot(roc[[28]], col="black", add=TRUE, lty =2)
plot(roc[[43]], col="red", add=TRUE, lty = 2)
q
knitr::opts_chunk$set(echo = TRUE)
# We will work within the Seurat framework
library(Seurat)
library(SeuratData)
data("pbmc3k")
pbmc <- pbmc3k
install.packages("reticulate")
install.packages("reticulate")
knitr::opts_chunk$set(echo = TRUE)
library(Seurat)
library(SeuratData)
library(reticulate)
reticulate::install_miniconda()
anndata::install_anndata()
install.packages("anndata")
library(anndata)
y=c(58.2,56.3,50.1,52.9,57.2,54.5,54.2,49.9,58.4,57.0,55.4,50.0,55.8,55.3,51.7,54.9)
catalyst=rep(c(1, 2, 3, 4), 5, 4, 3, 4) # matches y
catalyst
plot(y ~ catalyst, xlim=c(1, 5), ylab="catalyst concentration (y)", xlab="Mixture concentration", main="Catalyst experiment" )  # always plot the data first
data = c(176,6.5,.88,.42,
176,9.5,.88,.25,
190,9.0,1.00,.56,
176,8.9,.88,.23,
200,7.2,1.00,.23,
167,8.9,.83,.32,
188,8.0,.94,.37,
195,10.0,.98,.41,
176,8.0,.88,.33,
165,7.9,.84,.38,
158,6.9,.80,.27,
148,7.3,.74,.36,
149,5.2,.75,.21,
163,8.4,.81,.28,
170,7.2,.85,.34,
186,6.8,.94,.28,
146,7.3,.73,.30,
181,9.0,.90,.37,
149,6.4,.75,.46)
data = matrix(data,ncol=4,byrow=T)
bwt  = data[,1]
lwt  = data[,2]
dose = data[,3]
y    = data[,4]
n    = length(y)
out=lm(y~ bwt+lwt+dose)
summary(out)
plot(out)
for(i in 1:3){
plot(1:n,infl$coef[,i],pch=19,type="h")
lines(1:n,rep(0,n),lty=3,col=2)
}
plot(1:n,st.res,type="h")
lines(1:n,rep(0,n),lty=3,col=2)
print(data[3,])
par(mfrow=c(1,1))
library(Giotto)
library(devtools)  # if not installed: install.packages('devtools')
library(remotes)  # if not installed: install.packages('remotes')
remotes::install_github("RubD/Giotto")
remotes::install_github("RubD/Giotto@cless")
#!/usr/bin/env Rscript
source('~/Documents/R/seurat_functions.R')
#!/usr/bin/env Rscript
source('/Users/zhenwang/Downloads/CancerData\ /code/seurat_functions_public.R')
install.packages('MSigDB')
install.packages("msigdbr")
library('MSigDB')
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("msigdb")
library('msigdb')
#!/usr/bin/env Rscript
source('/Users/zhenwang/Downloads/CancerData\ /code/seurat_functions_public.R')
library(MSigDB)
devtools::install_github("mw201608/msigdb")
Result = msigdb.genesets(sets = c('c5.go.bp','c5.go.cc', 'c5.go.mf'), type = 'symbols', species = 'human')
#### Arguments ####
args = commandArgs(trailingOnly = TRUE)
if (length(args) == 3){
cancer = args[1]
species = args[2]
author = ''
sample = args[3]
} else if (length(args) == 4){
cancer = args[1]
species = args[2]
author = args[3]
sample = args[4]
}
numi.min = -Inf
ngene.min = -Inf
mt.max = Inf
range = 2:25
gmin = 5
# Combining
orig.ident = paste0(cancer, species, author, sample)
print(orig.ident)
#### Load ####
st = Load10X_Spatial(getwd())
library(Matrix)
library(data.table)
library(Seurat)
library(dplyr)
library(SPOTlight)
library(igraph)
library(RColorBrewer)
library(SeuratDisk)
library(future)
library(anndata)
# check the current active plan
plan("multiprocess", workers = 8)
options(future.globals.maxSize= 3*1024*1024^2)
start_time <- Sys.time()
# args<-commandArgs(T)
# scrna_path = args[1]
# spatial_path = args[2]
# celltype_final = args[3]
# output_path = args[4]
sc_file_path = '/Users/zhenwang/Documents/ST_analysis/Datasets/simulate_data_seqFISH_plus/output/scrna.h5ad'
spatial_file_path = '/Users/zhenwang/Documents/ST_analysis/Datasets/simulate_data_seqFISH_plus/output/spatial.h5ad'
output_file_path = '/Users/zhenwang/Documents/ST_analysis/Results/simulated_seqFISH_plus/'
Convert('/Users/zhenwang/Documents/ST_analysis/Datasets/simulate_data_seqFISH_plus/output/scrna.h5ad', dest = "h5seurat", overwrite = TRUE)
sc <- LoadH5Seurat(sc_file_path)
st <- LoadH5Seurat(spatial_file_path)
devtools::install_github('YingMa0107/CARD')
yes
devtools::install_github('YingMa0107/CARD')
devtools::install_github('YingMa0107/CARD')
install.packages('SingleCellExperiment')
version
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleCellExperiment")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleCellExperiment", force=True)
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleCellExperiment", force=TRUE)
version
library(SingleCellExperiment)
devtools::install_github('YingMa0107/CARD')
pkgbuild::check_build_tools(debug = TRUE
)
options(buildtools.check = function(action) TRUE )
R.version
version
install.packages('BiocManager')
install.packages('devtools')
devtools::install_github('YingMa0107/CARD')
BiocManager::install('SingleCellExperiment')
sessionInfo()
package.version('SingleCellExperiment')
BiocManager::install("SingleCellExperiment")
BiocManager::install("SingleCellExperiment", force=TRUE)
package.version('SingleCellExperiment')
sessionInfo()
library('SingleCellExperiment')
devtools::install_github('YingMa0107/CARD')
devtools::install_github('YingMa0107/CARD')
library(CARD)
install.packages('remotes')
install.packages("remotes")
library(devtools)  # if not installed: install.packages('devtools')
library(remotes)
remotes::install_github("RubD/Giotto")
library(Giotto)
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SPOTlight")
library(Giotto)
library(Giotto)
getwd()
my_python_path= "python path"
instrs = createGiottoInstructions(python_path = my_python_path)
my_python_path= "/Users/zhenwang/opt/anaconda3/bin/python3"
instrs = createGiottoInstructions(python_path = my_python_path)
sc_matrix<-read.table("../Datasets/simulate_data_seqFISH_plus/raw_somatosensory_sc_exp.txt",header=T,row.names = 1)
setwd("/Users/zhenwang/Documents/ST_analysis/Datasets/simulate_data_seqFISH_plus/")
my_python_path= "/Users/zhenwang/opt/anaconda3/bin/python3"
instrs = createGiottoInstructions(python_path = my_python_path)
sc_matrix<-read.table("./raw_somatosensory_sc_exp.txt",header=T,row.names = 1)
sc_lable<-read.table("./somatosensory_sc_labels.txt",header=F)
sc_matrix
print(dim(sc_matrix))
print(dim(sc_matrix))
print(dim(sc_label))
print(dim(sc_matrix))
print(dim(sc_lable))
sc_matrix[1,1]
sc_matrix[1,4]
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
install.packages("remotes")
remotes::install_github("satijalab/seurat-data")
library(SeuratData)
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
InstallData("stxBrain")
install.packages("http://seurat.nygenome.org/src/contrib/stxBrain.SeuratData_0.1.1.tar.gz", repos = NULL, type = "source")
install.packages("http://seurat.nygenome.org/src/contrib/stxBrain.SeuratData_0.1.1.tar.gz", repos = NULL, type = "source")
brain <- LoadData("stxBrain", type = "anterior1")
?stxBrain
library(COP)
library(Matrix) # isSymmetric
library(matrixcalc) # is.positive.definite
library(expm) # sqrtm
library(MASS)
library(dr)
library(rsample)
library(dplyr)
library(Giotto)
library(devtools)  # if not installed: install.packages('devtools')
library(remotes)  # if not installed: install.packages('remotes')
remotes::install_github("RubD/Giotto")
library(Giotto)
x = c(1993:2010)
y = c(2, 2, 1, 0, 0, 0, 0, -1, -1, 0, 0, 0, 1, 1, 1, 2, 2, 1)
plot(x,y)
plot(x, y, type = "l", col = "blue", lwd = 2, xlab = "X-axis", ylab = "Y-axis", main = "Line Plot Example")
points(node_x, node_y, col = "red", pch = 19)
plot(x, y, type = "l", col = "blue", lwd = 2, xlab = "X-axis", ylab = "Differece of number of communities", main = "Line Plot Example")
abline(h = horizontal_line_y, col = "red", lty = 2)
abline(h = 0, col = "red", lty = 2)
abline(h = 0, col = "blue", lty = 2)
plot(x, y, type = "l", col = "blue", lwd = 2, xlab = "Year", ylab = "After-Beofre", main = "Differece of number of communities")
abline(h = 0, col = "blue", lty = 2)
x = c(1994:2010)
y = c(0, 0, 1, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0)
plot(x, y, type = "l", col = "blue", lwd = 2, xlab = "Year", ylab = "After-Beofre", main = "Differece of number of communities")
abline(h = 0, col = "blue", lty = 2)
y = c( 0, 1, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0)
plot(x, y, type = "l", col = "blue", lwd = 2, xlab = "Year", ylab = "After-Beofre", main = "Differece of number of communities")
abline(h = 0, col = "blue", lty = 2)
x = c(1995:2010)
y = c( 1, 1, 1, 0, 0, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0)
plot(x, y, type = "l", col = "blue", lwd = 2, xlab = "Year", ylab = "After-Beofre", main = "Differece of number of communities")
abline(h = 0, col = "blue", lty = 2)
install.packages("pacman")
library(pacman)
pacman::p_load(phyloseq, metacoder, pryr, biomformat, RColorBrewer, ggplot2, gplots, Cairo, igraph, BiocParallel, randomForest, metagenomeSeq, MASS, DESeq2, vegan, RJSONIO, ggfortify, pheatmap, xtable, genefilter, data.table, reshape, stringr, ape, grid, gridExtra, splitstackshape, edgeR, globaltest, R.utils, viridis, ggrepel, ppcor, qs, dplyr, limma, memoise, tidyverse)
pacman::p_load(phyloseq, metacoder, pryr, biomformat, RColorBrewer, ggplot2, gplots, Cairo, igraph, BiocParallel, randomForest, metagenomeSeq, MASS, DESeq2, vegan, RJSONIO, ggfortify, pheatmap, xtable, genefilter, data.table, reshape, stringr, ape, grid, gridExtra, splitstackshape, edgeR, globaltest, R.utils, viridis, ggrepel, ppcor, qs, dplyr, limma, memoise, tidyverse)
pacman::p_load(phyloseq, metacoder, pryr, biomformat, RColorBrewer, ggplot2, gplots, Cairo, igraph, BiocParallel, randomForest, metagenomeSeq, MASS, DESeq2, vegan, RJSONIO, ggfortify, pheatmap, xtable, genefilter, data.table, reshape, stringr, ape, grid, gridExtra, splitstackshape, edgeR, globaltest, R.utils, viridis, ggrepel, ppcor, qs, dplyr, limma, memoise, tidyverse)
git clone https://github.com/xia-lab/MicrobiomeAnalystR.git
library(pacman)
pacman::p_load(phyloseq, metacoder, pryr, biomformat, RColorBrewer, ggplot2, gplots, Cairo, igraph, BiocParallel, randomForest, metagenomeSeq, MASS, DESeq2, vegan, RJSONIO, ggfortify, pheatmap, xtable, genefilter, data.table, reshape, stringr, ape, grid, gridExtra, splitstackshape, edgeR, globaltest, R.utils, viridis, ggrepel, ppcor, qs, dplyr, limma, memoise, tidyverse)
yes
install.packages("devtools")
library(devtools)
# Step 2: Install MicrobiomeAnalystR WITHOUT documentation
devtools::install_github("xia-lab/MicrobiomeAnalystR", build = TRUE, build_opts = c("--no-resave-data", "--no-manual", "--no-build-vignettes"))
devtools::install_github("xia-lab/MicrobiomeAnalystR", build = TRUE, build_opts = c("--no-resave-data", "--no-manual"))
vignette(package="MicrobiomeAnalystR")
library(MicrobiomeAnalystR)
sessioninfo()
version()
version
install.packages("pacman")
library(pacman)
pacman::p_load(phyloseq, metacoder, pryr, biomformat, RColorBrewer, ggplot2, gplots, Cairo, igraph, BiocParallel, randomForest, metagenomeSeq, MASS, DESeq2, vegan, RJSONIO, ggfortify, pheatmap, xtable, genefilter, data.table, reshape, stringr, ape, grid, gridExtra, splitstackshape, edgeR, globaltest, R.utils, viridis, ggrepel, ppcor, qs, dplyr, limma, memoise, tidyverse)
library(devtools)
devtools::install_github("xia-lab/MicrobiomeAnalystR", build = TRUE, build_opts = c("--no-resave-data", "--no-manual", "--no-build-vignettes"))
getwd()
data = read.txt("/Users/zhenwang/Downloads/real1.txt")
data = read.csv("/Users/zhenwang/Downloads/real1.txt")
data
dim(data)
library(devtools)
create("MultiCOP")
load_all("MultiCOP")
check("MultiCOP")
document("MultiCOP")
install("MultiCOP")
library(MultiCOP)
data = read.csv("/Users/zhenwang/Downloads/real1.txt")
data
getwd()
load("/Users/zhenwang/Library/CloudStorage/GoogleDrive-zwang9898@gmail.com/My\ Drive/project/dca_net/Data/dat_0609.RData")
dat
head(dat)
library(igraph)
library(sna)
library(reshape2)
library(HiveR)
library(SDMTools)
library(RColorBrewer)
library(plyr)
library(ggplot2)
library(GGally)
library(grid)
library(gridExtra)
library(reshape2)
library(igraph)
library(dplyr)
library(countrycode)
library(EnvCpt)
packageVersion("igraph")
packageVersion("reshape2")
packageVersion("ggplot2")
packageVersion("countrycode")
packageVersion("EnvCpt")
rm(list=ls())
getwd()
setwd("/Users/zhenwang/Documents/new_IR/example/")
knitr::opts_chunk$set(echo = TRUE)
# locate working dictionary on the example folder
source('../R/functions.R')
getwd("/Users/zhenwang/Documents/new_IR/TMCPD/example")
setwd("/Users/zhenwang/Documents/new_IR/TMCPD/example")
knitr::opts_chunk$set(echo = TRUE)
# set working dictionary on the example folder
source('../R/functions.R')
# set working dictionary on the example folder
source('../code/functions.R')
# Input: A sequence of networks represented in the form of adjacency matrices.
adj = load("adj_mat_seq.RData")
adj
load("adj_mat_seq.RData")
nonzero.dca
adj = nonzero.dca"
)
""
adj = nonzero.dca
save(adj, file="./adj_mat_seq.RData")
# set working dictionary on the example folder
source('../code/functions.R')
# Input: A sequence of networks represented in the form of adjacency matrices.
load("adj_mat_seq.RData")
yr = seq(1990,2010)
res = mod_cluster(Adj = nonzero.dca, yr, start=1)
library(EnvCpt)
fit_envcpt = envcpt(mod)  # Fit all models at once
mod = res$modularity
cl = res$membership
plot(c(1990:2010),mod,xlab="year",ylab="modularity")
library(EnvCpt)
fit_envcpt = envcpt(mod)  # Fit all models at once
fit_envcpt$summary  # Show log-likelihoods
plot(fit_envcpt)
fit_envcpt = envcpt(mod, models="meancpt")
fit_envcpt = envcpt(mod, models="meancpt")
plot(fit_envcpt)
cp.location = fit_envcpt$meancpt@cpts
cat(paste0("the 1st change point is", cp.location[1]))
library(EnvCpt)
fit_envcpt = envcpt(mod)  # Fit all models at once
fit_envcpt$summary  # Show log-likelihoods
fit_envcpt = envcpt(mod, models="meancpt")
plot(fit_envcpt)
cp.location = fit_envcpt$meancpt@cpts
cat(paste0("the 1st change point is at ", cp.location[1]))
?cat()
library(EnvCpt)
fit_envcpt = envcpt(mod)  # Fit all models at once
fit_envcpt$summary  # Show log-likelihoods
fit_envcpt = envcpt(mod, models="meancpt")
plot(fit_envcpt)
cp.location = fit_envcpt$meancpt@cpts
cat(paste0("the 1st change point is at ", cp.location[1]), "\n")
library(EnvCpt)
fit_envcpt = envcpt(mod, models="meancpt")
plot(fit_envcpt)
cp.location = fit_envcpt$meancpt@cpts
cat(paste0("the 1st change point is at ", cp.location[1]), "\n")
